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Abstract — Markov population models (MPMs) are a widely 
used modelling formalism in the area of computational biology 
and related areas. The semantics of a MPM is an infinite- 
state continuous-time Markov chain. In this paper, we use 
the established continuous stochastic logic (CSL) to express 
properties of Markov population models. This allows us to 
express important measures of biological systems, such as 
probabilistic reachability, survivability, oscillations, switching 
times between attractor regions, and various others. Because of 
the infinite state space, available analysis techniques only apply 
to a very restricted subset of CSL properties. We present a full 
algorithm for model checking CSL for MPMs, and provide 
experimental evidence showing that our method is effective. 

I. Introduction 

In the context of continuous-time Markov chains 
(CTMCs), properties of interest can be specified using con- 
tinuous stochastic logic (CSL) |[l], ||2l- CSL is a branching- 
time temporal logic inspired by CTL fS\. It allows to reason 
about properties of states (state labels), like the number of 
certain molecules given in this state, about what may happen 
in the next state (next operator), what may happen within 
a certain time (bounded until), what may finally happen 
(unbounded until) or about the long-run average behavior 
of a model (steady state). Because the underlying semantics 
of a Markov population model is given as a CTMC, we can 
also use CSL to reason about properties of such models, 
when interpreting CSL formulae on the CTMC semantics. 

We consider the complete set of CSL formulae, includ- 
ing the steady-state operator and in certain cases also the 
unbounded until operator fl\. The resulting logic can ex- 
press (nested) probabilistic properties such as "the long-run 
probability is at least 0.4 that we reach "^-states along $- 
states within time interval [6.5,8.5] with a probability larger 
than 0.98" via 5>o.4(^>o.98($Wf®-^'*-^' *))■ Using CSL, we 
can express many measures important for biological models, 
including oscillation I?)- 

Previous works Q, |l6l, Q have already considered 
techniques for the transient analysis of infinite-state CTMCs. 
These techniques are based on truncation. This means that 
only a finite relevant subset of the states of the infinite 
CTMC is taken into account. The extent to which these 
models are explored depends on the rates occurring there, 
as well as on the time bound of the transient analysis. Using 
truncation, techniques for the analysis of finite CTMCs [8] 
can be used for the analysis of properties of infinite-state 
models. 



In a previous publication ||9l, we have extended these 
results such that we were able to do approximate model 
checking for a subset of CSL. This subset excluded the 
steady-state operator as well as the unbounded until oper- 
ator. We have implemented these techniques in the model 
checker INFAMY [10|. On the other hand, we recently 
developed means to find subsets of states of Markov pop- 
ulation model CTMCs which contain the relevant steady- 
state probability mass [ 1 1 J . For each of these states, we 
also obtain lower and upper bounds of the steady state 
probability. These techniques have been implemented in the 
tool Geobound [12J. 

In this paper, we combine GEOBOUND and INFAMY, such 
that we can also handle the CSL steady-state operator. In 
addition, we introduce advanced truncation techniques which 
allow us to explore the model in a more advanced way, lead- 
ing to a smaller number of states being necessary to check 
properties. By also taking into account not only the time 
bound of CSL properties, but also the atomic propositions, 
we can further restrict the state space to be explored. In 
certain cases, this also allows us to handle the unbounded 
until operator. Using a ternary logic ifTsll . lfT4l . in contrast 
to previous publications, we can compute safe lower and 
upper bounds for probabilities. In turn, we can decide exactly 
whether a certain formula holds, does not hold or whether 
this cannot be decided on the finite truncation of the current 
model. Apart from this, we have also made some technical 
improvements, applying for instance to the portability and 
robustness of INFAMY. We show the applicability of the 
approach on a number of biological models. 

Organization of the paper: We give background on 
Markov population models, CTMCs and CSL in Section Ull 
and also recall the established CSL model checking algo- 
rithm for finite CTMCs. In Section Hill we give the main 
contribution of the paper, CSL model checking for infinite 
CTMCs. Section HVl reports experimental results. Section IVl 
gives related work and Section IVll concludes. 

II. Preliminaries 

A. Ternary Logic 

We consider a ternary logic ITSl . lfT4l with values := 
{T,±,?}. With the ordering _L < ? < T, B3 forms 
a complete lattice. We interpret n as the meet ("and" 
operator), and -"^ as the complement ("not") operation, with 
the usual definitions. Other operators Hke U ("or" operator) 
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Truth values of ternary logic operations. 



can be derived. Then, T and ± can be interpreted as values 
definitely true or false respectively, and ? is interpreted as an 
unknown value. We give an overview of the truth values in 
Table |T] Consider a formula over a number of values some 
of which are ?. If the value of this formula is different from 
?, we know that when inserting _L or T instead of some of 
these values, the result would still be the same. This way, 
in some cases we can obtain a truth value even though we 
do not known the truth values of some formula parts. 

B. Markov Chains 

We define some basics of Markov chains and CSL. Let 
AP denote a set of atomic propositions. For a countable set 
S, a distribution /i over 5 is a function /i : S* — > [0, 1] with 
fj.{S) = 1, and we define ^(A) :— J2seA ACS. 

Definition 1: A labelled continuous-time Markov chain 
(CTMC) is a tuple C = {S, Sinu, R, L) where 5 is a count- 
able set of states, Sinit is an initial state, L : {Sx AP) —> B3 
is a labelling function, and H : {S x S) ^ K>o is the rate 
matrix. 

In the stochastic process which the CTMC represents, Smit 
is the state we start in. The labelling function tells us for each 
state s and atomic proposition a whether a holds in s, does 
not hold in s, or whether this is not known or not specified. 
We say that the rate matrix is rate-bounded if the supre- 
mum supsg5 R(s, S) is finite, otherwise, it is called rate- 
unbounded. If R(s, s') > 0, we say that there is a transition 
from s to s'. For s e S, let post(s) := {s' \ R(s, s') > 0} 
denote the set of successor states of s. A state s is called 
absorbing if post(s) = 0. A CTMC is finitely branching if 
for each state s the set of successors post(s) is finite. In 
this paper, we consider rate-unbounded, finitely-branching 
CTMCs which do not explode l,15J . 1,161. Roughly speaking, 
if the CTMC explodes, then there is a positive, non-zero 
probability of infinitely many jumps in finite time. On the 
contrary, the fact that a CTMC does not explode implies that 
in finite time with probability one only a finite number of 
states can be reached. 

Transition probabilities in CTMCs are exponentially dis- 
tributed over time. The probability that an arbitrary transition 
is triggered within time t is given by 1 — e~^(*'''^)*, where 
R(s,yl) := Es'sA^-lS'*') f™" ACS. The probability of 
taking a particular transition from s to s' within time t is 
^'■'^''^ (1 — e^^'^^'^^^'Y State s' is reachable from s if there 



For a set of states A, by C[A] we denote the CTMC 
in which states of A have been made absorbing. For C = 
(S", Sj™t,R,i) and A C 5wehaveC[A] := (S*, s„,t, R', L) 
with R'(s,s') := R(s,s') for s ^ A and R'(s,s') := 
else, for all s,s' E S. Given a rate matrix R, we define the 
corresponding infinitesimal generator matrix Q such that 
for s, s' \f s ^ s' then Q(s, s') :— R(s,s') and we let 

Paths and probabilistic measure: Let 
C = (S*, Sinit, R, L) be a CTMC. A path is a sequence 
a = 51^152^2 • • • satisfying R(si, Si+i) > 0, and ti e R>o- 
Paths are either infinite or have a last state s„ with 
R(s„,5) = 0. For the path a and i G N, let a[i] — Si 
denote the i-th state, and let S{(T,i) — ti denote the time 
spent in Si. For t £ M>o, let a@t denote cr[i] such that i 
is the smallest index with t < Po'" PathP 

denote the set of all paths, and Path'''{s) denote the set of 
all paths starting from s. A probability measure Prf on 
measurable subsets of Path'''{s) is uniquely defined. We 
omit the superscript C if it is clear from context. 

C. Markov Population Models 

In this paper, we will consider Markov population models 
(MPMs), a sub-class of CTMC, where states encode the 
number of individuals of certain population types repre- 
sented by a vector over the natural numbers. Transitions 
between those states are defined by a change vector that 
characterizes the successor state and a rate determined by a 
propensity function that is evaluated in the predecessor state. 

Formally, a MPM with d population types is a CTMC 
C = (S',s™t,R,L) with S C N''. The rate matrix R will 
be induced by transition classes. 

Definition 2 (Transition Class): A transition class is a 
tuple (a, v) where a : N'' — > M>o is the propensity function 
and u e Z*^ \ {0} is the change vector. 
For the propensity functions, we will restrict to multivari- 
ate polynomials in N''. Given a set of transition classes 
{ti , . . . , Tm}, we define the entries of R via R(x, x + v) = 

E{j I v,=v}(^]i^) for '^j = ("i'^j) and I < j < m. 
In order to ensure R is well defined, we demand that for 
all states x £ S and all transition classes {aj,Vj) of our 



model aj{x) > only holds if 



e S. For the 



R(s,S) 

exists a sequence of states si 
Sn = s', and for each i = 1, . . 



. . ,Sn with n > 1, si = s, 
, n — 1 it is R(si, s^+i) > 0. 



labelling L we demand that each state is labelled by those 
boolean expressions over population counts and constants 
that evaluate to true. We also refer to the populations by their 
name. For example, if we have a MPM with populations A 
(xq) and B (xi), then state (3,4) would be labelled by e.g. 
A > c for all c > 3, ^ + B < 7 and + 52 <- 39 to name 

just a few. 

D. CSL Model Checking 

We consider the logic CSL jJJ interpreted over a ternary 
logic US, m. Let / = [t,t'] be an interval with t,t' C 
M>o U {00} with t' = 00 ^ t = and i < i'. Let p e [0, 1] 



Is, a} := Lis, a), [s.^iA^al := [s, -I'll □ [s, ^a], := 
[^,*]^ 

(t, PrsWlla, ^=T}c<p a Prs{(T\la, 4,}j^±}c<p, 
ls,V^pm:=l±, Prs{a\ia,4>l=T}:^p A PrsW\la,4,l^±}:^p, 
I ? else, 

|T, Sl(^)(^P a Su{^)txip, 
Is,5mp(*)1:=K, Sz,(*)0pA5c/('I>)^p, 
I ? else, 

where 

5l($) := lim Prs{a \ la@t,<S>} = T}, 

t—^OC 

Sui^) ■■= lim Prs{a | [aOt, $] ^ ±}, 

t — >OC 

h,^'"^} ■■= (T iff 5{a,0) el,± else) n [o-[l],*I, 

{T, BteLfam, $2l=T A Vi'<t.[[o-@t', *iI=T, 
±, vtei.lamt, <S>2}=± V 3t'<t.la@t', *il=±, 
? else. 

Table II 
CSL SEMANTICS. 



and M £ {<,<,>,>}■ The syntax of state formulae ($) 
and path formulae (</>) is: 

$ = a I I $ A $ I T'mpI'/') I S^p{<i>), 

Let F be the set of all CSL formulae. The truth value 
[•]c : ((5* U Path) x F) ^ B3 of formulae is defined 
inductively in Table If the model under consideration is 
clear from the context, we leave out the index C. The model 
C satisfies a formula $ if the initial state Sinit does. We 
specify the unbounded until as $1 U $2 := 'i'l $2 
and the eventually operator as := true $ where 

true = a V -la. Let i? C 5 be a set of states. For better 
readability, we use the name of the set B as an atomic 
proposition in formulae to characterize that the system is 
in a state contained in B. 

Model checking CSL formulae without the operators X 
and 5 on a finite CTMC over a ternary logic has already 
been used before to handle a different abstraction technique 
for CTMCs [131, flTl. Adding routines to handle X and S 
is straightforward. 

III. Model Checking based on Truncations 

Now we discuss how to model check CSL formulae on 
infinite CTMCs. To this end, our goal is to operate on a 
finite truncation instead of the original infinite CTMC. In 
a nutshell, starting from the initial state, we explore the 
states of the infinite model until we assume that we have 
explored enough of them to obtain a useful result in the 
following steps. Then, we remove all transitions from the 
states at the border and set all of their atomic propositions to 
?. In previous works |9| we already discussed some variants 
of such a model exploration. In Subsection IIII-AI we give 
another such technique, which is more efficient, because it 
needs to explore less states. Afterwards, we discuss how to 
build a finite CTMC truncation for a nested CSL formula. 



Finally, in Subsection llll-Cl we explain how to obtain results 
for the infinite-state CTMC only using the finite submodel. 

A. Truncation-based Reachability Analysis 

Given a set of states Wq, a CTMC C = (S*, s,™*, R, i) 
and CSL formula of the form $ ~ 'P^p{^i U $2), we want 
to compute a finite submodel C|yy = (WU W, R>v,L>v) 
sufficient to decide $ on all states of Wq. We define finite 
truncations of a CTMC. 

Definition 3: Let C = (5, s™t,R, L) be a CTMC. Let 
W C be a finite subset of S, and let W := post(yV) \ W. 
The finite truncation of C is the finite CTMC C\y^ := 
(WuW,Rw,iw) where Lw{s,-) := L{s,-) if s G W, 
and Lw(s,-) = ? else. The rate matrix is defined by 
Rw(s, s') := R(s, s') if s G W, and Rw(s, ■) := else. 

We build the truncation of the model iteratively, using 
the high-level (transition class) description of the model. 
Starting from states A, we explore the model until for all 
s £ Wo the probability to reach states in W is below an 
accuracy e, which we may choose as a fixed value or due 
to the probability bound p. 

Algorithm [T] describes how we can obtain a sufficiently 
large state set W. For s e S, s' eW and t e M>o U {00} 
we use 7Tc{s,t,s') to denote the probability that at time 
t G M>o, the CTMC C is in state s', under the condition 
that it was in s initially. For i = cx), we let tt denote 
the limit for < — !■ cx). As s' is absorbing, this value 
exists. Further, for a set of absorbing states B, we let 
S^c{s,t,B) := J2s'eB'^c{s,t,s') denote the probability to 
reach B within time t G R>o U {00}. Given a fixed s and 
t, we can compute 7rc(s, t, s') for all s' at once effectively, 
and given B and t we can compute £^{s,t,B) for all s at 
once effectively [2|. 

The algorithm is started on a CTMC C and a set of states 
Wq, for which we want to decide the property. We also 
provide the time bound t as well as the accuracy e. With W 
we denote a set of states for which the exploration algorithm 
may stop immediately, as further exploration is not needed 
to decide the given property. For $ above and / = [0, a], a G 
R>o U {00}, we could specify VV as the states which fulfill 
$2 V (-i^i A -'<i>2)- For all paths of the model, the truth 
value is independent of the states after such a state. 

B. Truncation-based Steady-state Analysis 

In the following, we will develop a technique to retrieve a 
finite subset of states W C S* that contains most of the total 
steady-state probability mass, i.e. J^cew ^('') > 1 ~ ^ for a 
given e < 1. The next step will be to derive lower and upper 
bounds on the state-wise steady-state probabilities inside that 
window W. 

Geometric Bounds: For the presented methodology to 
be applicable, we have to restrict our models to ergodic 
MPMs, since for steady-state analysis the equilibrium dis- 
tribution has to exist uniquely. Ergodicity can be verified by 



Algorithm 1: TransientTrunc(C, Wo, VV, e)- 
W Wo 

W := post(W) \ W 

while maxsewo ^c[w] ^, W \ W) > e do 

choose s from arg maxsg yVo fc[W] ^ \ ^) 
while ^c[W] i, W \ W) > e do 
choose A C (W \ W) such that 

W := WUA 
W := post(W) \ W 
return W U (post(W) n W) 



the means of Lyapunov functions and the following theorem. 
In the following, by X{t) we refer to the stochastic process 
underlying the MPM, and by E to the expectation of a 
random variable. 

Definition 4 (Lyapunov Function): A Lyapunov function 
is a function 5 : N'' — M>o. 

Theorem 1 (Tweedie (JTj): Assuming that Xit) is irre- 
ducible, it is also ergodic and uniquely determined by its 
infinitesimal generator iff there exists a Lyapunov function 
g*, a finite subset W C 5, and a constant A > such that 

1) ^E[g*{X{t)) I X{t) = x]<-\ for all a: e N''\ W, 

2) \E[g*{X{t)) I X{t) x] < 00 for all a; e W, 

3) the set {x e N"* | g*{x) < 1} is finite for all / < 00. 
Therefore, given a MPM C = {S, Si„if , R, L) induced by a 
set of transition classes {ti, . . . ,t„i}, we can use Theorem 
[1] for a semi-decision procedure to check ergodicity by 
choosing candidates for Lyapunov functions. Our experience 
has shown that in most cases rather simple functions, i.e. 
multivariate polynomials of degree two, already suffice for 
usual models from systems biology and queuing theory. 
Consequently, we restrict the choice of Lyapunov functions 
to that class. 

We can exploit Theorem [T] also to retrieve the aforemen- 
tioned window W that encloses most of the steady-state 
probability mass. 

At first, let the drift d* (x) in state x e N'' be defined as 

d*{x) := ±E[g*{X(t)) \ X{t) ^ x] ^ W){x). 

Due to the transition class induced structure of Q, the drift 
is given as 

m 

d*{x) ^^aj{x){g*{x + Vj) -g*{x)) 

and can easily be represented symbolically in a state x. 

Since the propensity functions of our models as well as 
the Lyapunov function are multivariate polynomials, so is 
the drift. The next step is to retrieve a positive real number 
c > maxj-gfijd In order to find that global maximal 



drift, common global optimization techniques like gradient 
based methods and simulated annealing can not be used, 
since there is no guarantee to get the real global maximum. 
What we propose instead is to solve '\Id*{x)\p = for all 
p £ - where f{x)\p denotes the projection of the 
multivariate function / onto the subspace spanned by a;^ e p 
to retrieve all K possible candidates rrifc. In order to solve 
the equation systems Vd*{x)\p = 0, we suggest the use 
of the polyhedral homotopy continuation method, which is 
guaranteed to find all roots. For implementations, we refer 
to ifTsl and lfT9l . Finally, after restricting the candidate set 
to M = {ruk I TOfc G K>o}, we set c = ma.XmeM d* {m). 
Please note, that due to The existence of a maximal value 
c, the chosen Lyapunov function serves as a witness for 
ergodicity assuming an irreducible MPM. By scaling the 
Lyapunov function by we retrieve the normalized drift 

dix) - ^^E[g{X{t)) I X{t) - x] 

and using conditions 1 and 2 from Theorem [T] we get 

c _ 

d{x) < — ■ Xw{x), (1) 

c + 7 

where xw(x) — 1 if c ^ W and else. Multiplying 
Inequality [U with 7t{x) and summing over x leads us to 

^(W)= ^4x)<-^. 

Convergence of this sum is guaranteed for (infinite) ergodic 
MPMs as stated in |20|. Consequently, we can exploit this 
inequality by directly choosing d{x) = ^d*{x), i.e. setting 
A > such that e = to get 7r(W) > 1 - £ for 

W = {x e N"* I d{x) >£-!}. 

State-wise Bounds: Given our state space window W, 
our next goal is to get lower {l{x)) and upper (u{x)) bounds 
on the steady-state probabilities inside W, i.e. probability 
vectors / and u such that l{x) < Tr{x) < u{x) for all c £ 
W. For this, we will employ the methodology developed by 
Courtois and Semal 1211 ll22l which we have extended to 
infinite state MPMs in ifTTI : 

Theorem 2 (RTTi): Let Q be the infinitesimal generator of 
an ergodic CTMC X{t) with countably infinite state space S 
and let W C S* be a finite subset of the state space. Further, 
we let the matrix C be the finite submatrix of Q containing 
exactly the states in W. If by U we refer to the uniformized 
CTMC of C, i.e. 

U = I + a^^C with a > max — C(«, z), 

i 

then for all x £ W we have 

min vr^j (x) < ^^^^ < max tt^J (x) 
1 Lxew ^(3^) ^ 



where tt j is the steady-state distribution of matrix Uj 
which is matrix U made stochastic by increasing column 
j- 

As presented in the previous paragraph, we have 1 — e < 
X^xew ''^{^) ^ 1- Consequently, we can use the geometric 
bounding technique to obtain the unconditional state-wise 
bounds from Theorem |2] that is to retrieve for all states 

X e w 

l{x) = (1 — e) min tt^J [x) < 7r{x) < maxTr^j (x) = u{x). 

j j 

Geobound: All these presented techniques, i.e., the 
retrieval of geometric bound via Lyapunov functions and 
the computation of state-wise steady-state bounds via the 
methodology of the previous paragraph have been imple- 
mented in a prototypical tool called GEOBOUND Iil2i . 

C. Truncation-based CSL Model Checking 

Given a CTMC C, we want to check whether C satisfies 
$. This is done in two phases. At first, we construct a 
finite truncation that is sufficient to check the formula. To 
this end, we employ an algorithm to determine suitable 
truncations. The states explored depend on the specific 
CSL formula analyzed. The computation works by recursive 
descent into sub-formulae. The most intricate formulae are 
the probabilistic operators, for which we use the technique 
from Section IIII-AI After the exploration, we can compute 
\sinit^ $1 on the finite truncation. 



Algorithm 2: Truncate(C = (5, Sj^t, R, L), $). CSL 
state space exploration. 

function Trunc(C, W, $) 
switch $ do 

case a return W 

case return Trunc(C, W, 

case $1 A $2 return 

Trunc(C, W, $i) U Trunc(C, W, $2) 

case V^piX^-^) return Trunc(C, WU W, 

case $2) 

Wt = 

TransientTrunc(C, W, stop(<I)), t, e) 

Wt' = 

TransientTrunc(C, Wf , stop($), t'-t, e) 
return 

Trunc(C, Wf, $1) U Trunc(C, Wt-, $2) 
case S^p{^) return 

TRUNC(C,pOSt(W)\ W,*) 
return C\^R^^c{c,{s,„t},'f) 



Algorithm|2]describes the exploration component. Given a 
CTMC C = (5, s™t, R, L) be a CTMC and a state formula 
we call Trunc(C, {sinit}, $). Afterwards, we can use 
the CSL model checking algorithm for a ternary logic on 
the model obtained this way. With stop($) we denote a set 



of states for which we can stop the exploration immediately, 
as exemplified in Section IIII-AI For nested formulae, this 
value is computed by a simple precomputation in a recursive 
manner. 

We employ ternary CSL model checking on the finite 
model obtained. However, we have akeady obtained the 
steady-state probabilities beforehand. Thus, to obtain the 
lower bound probabilities of 5|xp(^'), we sum up the lower 
bound steady-state probabilities of states s £ W with 
[x, = T. For the upper probability bound, we sum upper 
steady-state probabilities of states s with '5] ^ -L and 
add the probability e that limits the steady-state probability 
outside W. The probabilities computed are the probabilities 
for all states, because the model is ergodic. 

Correctness: Consider a truncation C|yy = 
(WuW, R>v,iw) constructed for the input CTMC 
C and a state formula $. If we obtain the truth value 
|s, $]] 7^? in Cyv, then this is also the truth value in C: The 
correctness is independent of the exploration algorithm, 
which plays a role for performance and applicability of 
the approach. If too many states are explored, we may run 
out of memory, but if too few are explored, we are unable 
to decide the value in the original model. As a result, the 
correctness of the algorithm for CSL without steady state 
follows by giving a simulation relation lfT4l Definition 
3.4.2] between C and and |14, Theorem 4.5.2]. The 
correctness of the steady-state extension follows as we give 
safe upper and lower bounds in exactly the same way as it 
is done in llT4l Theorem 4.5.2]. 

IV. Experimental Results 

Using several case studies, we assess the effectiveness of 
our technique. We have combined the GEOBOUND tool lfT2l 
to compute bounds on steady-state probabilities for Markov 
population models with the infinite-state model checker 
Infamy ifTOl . This way, we can effectively handle the 
combination of models and properties described in this 
paper. To show the efficiency of the approach, we applied 
our tool chain on a number of models from different areas. 
The results were obtained on an Ubuntu 10.04 machine with 
an Intel dual-core processor at 2.66 GHz equipped with 
3 GB of RAM. The tools used are available'. Instead of 
the truth value for the formula under considerations, in the 
result tables we give intervals of the probability measure of 
the outer formula. We also make use of a derived operator 
for conditional steady-state measures defined as 

{T, 5L(<I>i|<I>2)[X]pAS't/(<I>i|<I>2)cXp, 
^, 5l($i|$2)i?^pA5c/($i|$2)^P, 
? else, 

' ihttp : / /alma ■ cs ■ uni-saarland. de/?page_id=7 4| 



where 



$2) 
$2) 



SlI'^i a $2) 
A$2) 



5l($2) 

Protein Synthesis l!23\l : We analyze the MPM encoding 
protein synthesis, as depicted in Table Hill In biological cells, 
each protein (P, X2) is encoded by a gene (G, xi). If the 
gene is active (G = 1), the corresponding protein will be 
synthesized with rate = 1.0. Proteins may degenerate with 
rate 6 = 0.02 and thus disappear after a time. The gene 
switches from active state to inactive (G = 0) with rate /i = 
5.0 and vice versa with rate A = 1.0. Note that in a previous 
paper, this model has been presented as a stochastic Petri 
net (SPN). Often, transition class models and SPN (without 
zero-arcs) can trivially be encoded within each other. 
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Table IE 

Transition classes of the protein synthesis model. 



We consider the property that on the long run, given that 
there are more than 20 proteins, a state with 20 or less 
proteins is most likely (with a probability of at least 0.9) 
reached withing t time units: 



5>p(P>o.9(0[°'*lP < 20) 



P > 20). 



We give results in Table |IV] The shortcut stop of Al- 
gorithm |2] was not used for the analysis. With "depth" we 
specify the number of states of the shortest path from the 
initial state to any other state of the truncation. The runtimes 
of Geobound and Infamy is given in seconds. The rate 
of decay of proteins depends on the number of proteins 
existing. For the states on the border of W, we have large 
rates back to existing states. Because of this, for the given 
parameters the state space exploration algorithm does not 
need to explore further states, and the total number states 
explored n does not increase with the time bound. To obtain 
different n, we would have needed to choose extremely 
large time bounds, for which analysis would be on the one 
hand infeasible and on the other hand would lead to results 
extremely close to 1. The lower and upper bounds are further 
away than e. This results, because we have to divide by 
811(^2) for the lower and by Sl{^2) for the lower bound. 
In turn, this may lead to a much larger error than e. 

Gene Expression l[24]l : Next, we analyze a network of 
chemical reactions where a gene is transcribed into mRNA 
(M) with rate A = 25.0 and the mRNA is translated into 



proteins (P) with rate fi = 1.0. Both populations can 
degrade at rates Sm = 2.0 and 6p = 1.0, respectively. The 
corresponding transition classes are listed in Table [V] The 
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Table V 

Transition classes of the gene expression model. 



property of interest is the steady-state probability of leaving 
a certain set of states W enclosing more than 80% of the 
steady-state probability mass most likely within t time units, 
i.e., 

'5>p(7'>o.9(0[°'*l-W^) I W), 



where 



:= M > 5 A Af < 20 A P > 5 A P < 20 



with 5>o.8(W^) = T. 

The results are stated in Table |VT] Similar to the protein 
synthesis case study, we see that there is no increase in the 
number of states, because the window size already comprises 
enough states for the transient analysis. In Table IVIII we 
consider results for the subformula P>o.9(0'°'*'~'W^). We 
compare the methods to explore states for the transient 
until described in this paper (Advanced) with the finite state 
projection Q (FSP) previously used in INFAMY. We see 
that the time needed is comparable, but the new algorithm 
needs to explore less states. This is the case because with the 
method introduced here when building the finite truncation 
we have more control into which direction we explore. In 
contrast, the FSP explores the model into all directions at the 
same time, until enough precision is reached. When we use 
the shortcut stop (Advanced-nAP), we will not explore states 
further in which -^W holds. When exploring the model, 
for larger time bounds there is some point at which there 
are almost only states of which all successors have been 
completely explored and states for which -^W holds. Thus, 
the maximal number of states to be explored is constant with 
this optimization, except for very large time bounds. 

For the protein synthesis, there is almost no difference in 
the number of states needed by the new method and FSP. 
Because it has only one infinite variable, there is just one 
direction to be explored. Thus, the new method performs 
worse, as it needs more effort to explore into this direction. 

Using the shortcut method also allows us to handle 
the formula VyQ.g{()^W), involving the unbounded until 
operator. Using a precision of 10^®, we needed a total 
time of 1.5 seconds, reached 974 states and obtained a 
reachability probability almost 1.0. Using a precision of 



Table IV 
Protein synthesis results. 
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Table VI 
Gene expression results. 
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Table VII 

Gene expression - comparison of methods. 
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IQ- , we needed 3.7 seconds and explored 1823 states. 
For the computation of unbounded until probabilities, effi- 
cient specialized algorithms were used, which explains that 
we needed less time than for some of the time bounded 
experiments. 

Exclusive Switch f25^: The exclusive switch is a gene 
regulatory network with one promotor region shared by 
two genes. That promotor region can either be unbound 
(G = l,G.Pi = 0,G.P2 = 0) or bound by a protein 
expressed by gene 1 (G ^ 0, G.Pi = 1,G.P2 = 0) or 
gene 2 (G = 0, G.Pi = 0, G.P2 = 1). If the promotor 
is unbound, both proteins are expressed, each with rate 
p = 0.05, otherwise only the protein that is currently bound 



to the promotor is produced at rate p. Proteins degrade at rate 
S — 0.005. Binding happens at rate A = 0.01 and unbinding 
at rate p — 0.008. The transition class structure is given in 
Table IVlIll 



This system has two attractor regions, i.e., two spatial 
maxima in the steady-state probability distribution over the 
protein levels, one at Pi = 10, P2 = and the other one at 
Pi = 0, P2 = 10. We are interested in the switching time 
between these two regions. For this, we estimate the time 
needed for a 90%-quantile of the steady-state probability 
mass of one of the two attractors to reach the other attractor 
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Table VIII 

Transition classes of the exclusive switch. 



region. More precisely, let 

start := | |(Pi, P2) - (10, 0)| |^ < 4 
end ■= ||(Pi,P2)-(0,10)||2 <4. 

Then, the formula to check is 

5>p(7'>o.9(0["'*lend) I start). 

Note that since the model is symmetric we only have to 
check one formula from one attractor to the other. The 
corresponding results are depicted in Table |IX] 

From these results we may conclude that in half of the 
cases, most likely the switching time between the attractor 
regions is at most 7800 time units, while in almost all cases 
the switching time is most likely below 8000 time units, 
assuming the system has stabilized to a steady state. 

V. Related Work 

The techniques of this paper are derived from combina- 
tions of our previous works fTO^, f9l and {TT\, |T2|. This 
work has been inspired and is related by a number of other 
works. 

Finite state projection (FSP) by Munsky and Kham- 
mash |7| is closely related. The method also works on 
building a finite truncation of the original model. The proofs 
given work for general truncations, but in their publications 
they always use an algorithm which explores the model in a 
breadth-first way. They consider time-bounded reachability 
and no logic like CSL. Adaptive uniformization for CTMCs 
was introduced by van Moorsel and Sanders |6|. In this 
approach uniformization is recalibrated to perform well 
when exploring the state space on the fly. Remke et al. ^261 
have developed algorithms for model checking CSL against 
infinite-state CTMCs of QBDs and JQNs. The systems to 
which the method is applicable are less general, but the 
approach is less expensive than our method. To the best of 
the authors knowledge, the method to bound the steady-state 
probabilities using Lyapunov functions is used here for the 
first time for general MPMs. 



VI. Conclusion 

In this paper, we have shown how to model check CSL 
on Markov population models of infinite size. Without the 
steady-state operator, the method is also applicable for gen- 
eral CTMC derived from a high-level specification. We have 
evaluated our method on models from the biological domain. 
The method extends previous related publications by means 
to check the steady-state operator and gives guarantees for 
the truth value obtained. 
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